#!/usr/bin/python
#coding=utf-8

import numpy
import random as random_number


#самый простой алгоритм Монте-Карло
def MCint(f,n,M,a,b): 
    V=(b-a).prod()
    fff=map(f,[[numpy.random.uniform(low=a[i], high=b[i]) for i in xrange(n)] for k in xrange(M)])
    In=V*reduce(lambda x,y: x+y , fff)/M
    Di=V*reduce(lambda x,y: x+y , map(lambda x:x*x,fff))/M-In**2

    return (In,Di)












#разбиение области - генерация сетки для Монте-Карло
def rebin(rc,nd,r,xin,xi):
    k=0   
    dr=0.0
    xn=0.0
    xo=0.0
    for i in xrange(nd):
        while rc > dr:
              
              dr += r[k]
              k+=1

        if k >1:
           xo = xi[k-1]
        xn=xi[k]
        dr -= rc
        xin[i]=xn-(xn-xo)*dr/r[k]

    xi[nd] = 1.0
#алгоритм вегаса:
#regn - область инт-ия, itmx - число итераций, ncall - число вызовов на каждой итерации, init - начальное значение интеграла, ndim - размерность пространства, fxn - интегранд, tgral - выходное значение: значение интеграла, sd - погрешность, chi2a - значение хи-квадрата.
def MC_vegas(regn,itmx,ncall,ndim,fxn):
    xi=numpy.zeros( (ndim,itmx+1) )
    #константы, отвечающие за адаптивность
    ALPH=1.5
    TINY=1.0e-12   


    for j in xrange(ndim):
          xi[j,1]=1.0
   
    si=swgt=schi=0.0
    

       

    
    

    dv2g = 1.0/(ncall-1.0)   
    xnd = 1.0
    xjac = 1.0/ncall
    dx=numpy.zeros(ndim)
    for j in xrange(ndim):
        dx[j]=regn[j+ndim]-regn[j]
        xjac *= dx[j]
    r=numpy.zeros(itmx+1)
    xin=numpy.zeros(itmx+1)
    
    for i in xrange(itmx):
        r[i]=1.0
    for j in xrange(ndim):
        rebin(1.0/xnd,itmx,r,xin,xi[j])

    kg = numpy.zeros(ndim)
    x=numpy.zeros(ndim)
    ia = numpy.zeros(ndim)
    d=numpy.zeros( (itmx+1,ndim) )
    di=numpy.zeros( (itmx+1,ndim) )
    dt = numpy.zeros(ndim)
    for it in xrange(itmx):
        ti=tsi=0.0
        for j in xrange(ndim):
            kg[j]=1.0
            for i in xrange(itmx):
                d[i,j]=di[i,j]=0.0
   
            
        fb=f2b=0.0
        for k in xrange(ncall):
            wgt = xjac
            for j in range(ndim):
                xn=(kg[j]-random_number.uniform(0,1))+1.0
                ia[j] = numpy.floor(max(min(xn,itmx),1))                      
                if ia[j] > 1:
                   xo = xi[j,ia[j]]-xi[j,ia[j]-1]
                   rc = xi[j,ia[j]-1]+(xn-ia[j])*xo
                else:
                   xo = xi[j,ia[j]]
                   rc = (xn-ia[j])*xo
                         
                   x[j]=regn[j]+rc*dx[j]
                  
                   wgt *= xo*xnd
                   
            f=wgt*fxn(x)
            if fxn(x) >100:
               print fxn(x), x        
            f2=f*f
            fb += f
            f2b += f2
                   
        for j in xrange(ndim):
            di[ia[j],j]+=f
            d[ia[j],j]+=f2

        f2b = numpy.sqrt(f2b*ncall)
        f2b = (f2b-fb)*(f2b+fb)
        ti += fb
        tsi += f2b
                 
                  
   
        tsi *= dv2g
        wgt=1.0/tsi
        si += wgt*ti
        schi += wgt*ti*ti
        swgt += wgt
        tgral=si/swgt
#        chi2a=(schi-si*tgral)/(it-0.9999)
#        if chi2a < 0.0:
#           chi2a = 0.0
        sd=numpy.sqrt(1.0/swgt)
        tsi=numpy.sqrt(tsi);
       

        for j in xrange(ndim):

            xo=d[1,j]

            xn=d[2,j]

            d[1,j]=(xo+xn)/2.0

            dt[j]=d[1,j]

            for i in range(2,itmx-1):
                rc=xo+xn
                xo=xn
                xn=d[i+1,j]
                d[i,j] = (rc+xn)/3.0
                dt[j] += d[i,j]

                d[itmx,j]=(xo+xn)/2.0
                dt[j] += d[itmx,j]

             


         
    return (tgral,sd)


